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ABSTRACT 

We present results of a pulsar population synthesis study that incorporates 
a number of recent developments and some significant improvements over our 
previous study. We have included the results of the Parkes multi-beam pulsar 
survey in our select group of nine radio surveys, doubling our sample of radio 
pulsars. More realistic geometries for the radio and 7-ray beams are included in 
our Monte Carlo computer code that simulates the characteristics of the Galactic 
population of radio and 7-ray pulsars. We adopted with some modifications the 
radio beam geometry of Arzoumanian, Chernoff & Cordes (2002). For the 7-ray 
beam, we have assumed the slot gap geometry described in the work of Muslimov 
& Harding (2003). To account for the shape of the distribution of radio pulsars 
in the P — P diagram, we continue to find that decay of the magnetic field on 
a timescale of 2.8 Myr is needed. With all nine surveys, our model predicts 
that EGRET should have seen 7 radio-quiet (below the sensitivity of these radio 
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surveys) and 19 radio-loud 7-ray pulsars. AGILE (nominal sensitivity map) is 
expected to detect 13 radio-quiet and 37 radio-loud 7-ray pulsars, while GLAST, 
with greater sensitivity is expected to detect 276 radio-quiet and 344 radio-loud 
7-ray pulsars. When the Parkes multi-beam pulsar survey is excluded, the ratio 
of radio-loud to radio-quiet 7-ray pulsars decreases, especially for GLAST. The 
decrease for EGRET is 45%, implying that some fraction of EGRET unidentified 
sources are radio-loud 7-ray pulsars. In the radio geometry adopted, short period 
pulsars are core dominated. Unlike the EGRET 7-ray pulsars, our model predicts 
that when two 7-ray peaks appear in the pulse profile, a dominant radio core 
peak appears in between the 7-ray peaks. Our findings suggest that further 
improvements are required in describing both the radio and 7-ray geometries. 

Subject headings: radiation mechanisms: non-thermal — magnetic fields — stars: 
neutron — pulsars: general — 7 rays: theory 


1. Introduction 

Rotation-powered pulsars are the brightest class of 7-ray sources detected by the Comp- 
ton 7- Ray Observatory (CGRO). The high-energy telescope EGRET made firm detections 
of pulsed 7-ray emission from five known radio pulsars (Thompson 2001), and possible de- 
tections from several others (Kanbach 2002). In addition, the high-energy pulsar Geminga 
may be radio-quiet, or at least radio-weak. EGRET also discovered more than 200 7-ray 
sources (Hartman et al. 1999), most of which are still unidentified. However, several dozen 
new radio pulsars, out of more than 600 discovered since the end of the CGRO mission by 
the recent Parkes multi-beam pulsar survey (PMBPS) (Manchester et al. 2001) or in deep 
targeted observations (Lorimer 2003), lie within the error circles of EGRET sources. Al- 
though many of these are young, energetic pulsars, their identification as 7-ray pulsars must 
await observation with the next 7-ray telescopes, AGILE and GLAST. 

In the meantime, population synthesis studies of radio and 7- ray pulsars can predict the 
number of radio-loud and radio-quiet 7-ray pulsar detections expected by different telescopes, 
assuming different models for radio and 7-ray emission. Even though current radio and 7-ray 
emission models have a number of outstanding uncertainties, the results of such studies can 
provide quite sensitive discrimination between models. In particular, polar cap and outer 
gap models of 7-ray emission make very different predictions of the number of radio-loud 
and radio-quiet 7-ray pulsars. Polar cap models (e.g. Daugherty & Harding 1996), where 
the high-energy emission region is located on the same open field lines as the radio emission, 
expect a large overlap in the radio and 7-ray emission beams and thus a higher ratio of radio- 
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loud to radio-quiet 7-ray pulsars. On the other hand, outer gap models predict a smaller 
overlap between 7-ray and radio emission beams, because the high-energy and visible radio 
emission originate from opposite poles (Romani & Yadigaroglu 1995, Cheng et al. 2000), 
thus predicting more radio-quiet than radio-loud 7-ray pulsars. Thus, population synthesis 
can also address the question of how many EGRET unidentified sources are radio pulsars. 

Results of our initial study of pulsars in the Galaxy were presented by Gonthier et al. 
(2002). In this work, we evolved neutron stars from birth distributions in space, magnetic 
field strength, period and kick velocity, in the Galactic potential to simulate the population 
of radio pulsars detected in eight surveys of the Princeton catalog (Taylor, Manchester & 
Lyne 1993). A very simple model of radio and 7-ray beams was assumed, in which both 
were aligned with solid angle of 1 sr. Radio luminosity was assigned using the model of 
Narayan & Ostriker (1990) and 7-ray luminosity from the polar cap model of Zhang & 
Harding (2000). We found that agreement of the distribution of simulated radio pulsars 
with the observed distribution was significantly improved by assuming decay of the neutron 
star surface magnetic field on a time scale of 5 Myr. With these assumptions, EGRET should 
have detected 9 radio-loud and 2 radio-quiet 7-ray pulsars, and GLAST should detect 90 
radio-loud and 101 radio-quiet pulsars (9 detected as pulsed sources). Because the radio 
and 7-ray beam apertures were assumed to be identical, radio-quiet 7-ray pulsars were those 
whose radio emission is too weak to be detected by the selected radio surveys. 

There have been a number of new developments in both radio pulsar observation and 
analysis, and in 7-ray pulsar theory since we completed our initial population study. The 
PMBPS (Manchester et al. 2001, Morris et al. 2002, Kramer et al. 2003) is nearly complete 
and has more than doubled the number of radio pulsars with measured period derivatives 
from 445 to nearly 1300. Determination of pulsar distances from dispersion measure has 
been greatly improved with the development of a new model by Cordes & Lazio (2002). New 
radio luminosity and beam models have been developed by Arzoumanian, Chernoff & Cordes 
(2002), which describe core and conal components of the emission and their dependence on 
period and frequency. Arzoumanian, Chernoff & Cordes (2002) have also derived a new 
two-component distribution of radio pulsar space velocities. A new polar cap 7-ray emission 
model has been developed by Muslimov & Harding (2003), in which radiation from pair 
cascades at high altitude along the edge of a slot gap forms a wide hollow cone of emission. 
In addition, the solid angle as well as the luminosity of the 7-ray beam is described in this 
model. 

This paper presents results of an expanded and updated pulsar population synthesis 
study that includes all of the above recent developments, as well as improved 7-ray sensitivity 
maps. By incorporating independent models for the radio and 7-ray beam geometry, we are 
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now able to investigate how the beam geometry affects the observable characteristics of radio- 
loud and radio-quiet 7-ray pulsar populations. We are also able to address the question of 
how many EGRET unidentified sources are expected to be radio-loud or radio-quiet 7-ray 
pulsars in the polar cap model. Of particular interest is the issue of how many of the new 
Parkes radio pulsars in EGRET error circles are counterpart 7-ray pulsars. In addition, we 
can make more accurate estimates of the numbers of radio-loud or radio-quiet 7-ray pulsars 
detectable by the AGILE and GLAST telescopes. 


2. Radio Emission Geometry 

We have adopted the geometry model for the radio emission beams as presented by 
Arzoumanian, Chernoff k Cordes (2002) (ACC from now on) with some slight modifica- 
tions. We have assumed a core and a single conal component described by Gaussians with 
characteristic widths as follows 


p core = 1.5 °P-°- 5 , and 

Pcone = 5.2°(l + ^P-°- 5 (1) 

where the period, P, is in seconds and the frequency, u, is in MHz. The characteristic core 
width is the width of the core beam at 1/e of the peak intensity. We have incorporated the 
radius-to-frequency mapping in the conal width developed by Mitra k Deshpande (1999); 
although they introduce elliptical shapes to the conal geometry, they find no compelling 
reason to abandon circular beams. The coefficient of 5.2° above is chosen to give the same 
width at 400 MHz as in ACC. 

For each simulated pulsar, the pulse profile is binned into 500 bins of the phase angle, 
(j>, ranging from -7r to 7 r. Each bin is assigned a flux, s((j), is), consisting of the sum of the 
flux contributions from the core and cone components given by 


s((j), v) 

Si(9,u) 


Score(9) v) H“ Sconei®) v) 


Oii + 1 


( •'-) 

v \50 MHz/ 


a=i +1 Jj 4 

mf, 


( 2 ) 


where i indicates core or cone, cc* is the spectral index (Si oc v ai with on < 0), 6 is the 
polar angle to the magnetic axis, /,(0) is the angular distribution of the component flux, L, 
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is the component luminosity and d is the distance to the pulsar. The relationship between 
the phase angle, (f>, and the polar angle, 9, depends on the viewing geometry given by the 
expression 


cos# = sin a sin C cos + cos a cos C- (3) 

During the simulation, the magnetic inclination angle, a, and the observer’s line of sight 
angle, £, are chosen randomly between zero and 7t/2, accounting for emission from both 
poles. The difference between these two angles defines the impact angle, f3 = £ — a. 

We assume that the spectrum of each component has a low frequency cutoff of 50 
MHz and can be modeled by a single power law, with spectral indices of a core = — 2.1 and 
otcone = — 1.6. It is general agreed that the spectra of cores are steeper that those of cones 
especially for short-period pulsars (Lyne k Manchester 1988 and Rankin 1983). As discussed 
later in section 5, we have assumed constant spectral indices with a difference of 0.5 between 
the core and cone indices. 

The angular distributions of the core and conal components are given by the Gaussians 


fcore(0) = -^—e and (4) 

^ ^core 

/«».(») = 

^ cone 

The solid angles for each of the components are chosen to normalize the Gaussian distribu- 
tions describing the angular distribution of the flux in equation (5) when integrating over 
the polar angle, 6, the azimuthal angle, <f>, and are given by the approximate expressions 


flcore — ^ P corei S' 11 *! (5) 

flcone = 2n 3/2 w e d ~ 0.&np 2 cone . 

The width and radius of the annulus of the conal beam are given by 


w e 
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These expression differ slightly from the ones used by Arzoumanian (private communication) 
given by the following forms: 


n 


ACC _ 
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ACC _ 


^Pcore 

In 2 

^][Pcont 

5 In 2 


1.47rp; 


, 2 , 

core’ 


and 


0.587rp; 


, 2 

cone' 


(7) 


With the above definition of w e when 9 = p cone , p cone represents the radius of the cone at 
half power. The half-angle coefficient of 3.9° in equation (6) is the angle where the conal 
intensity peaks. This coefficient corresponds to that of the middle cone of the three conal 
components discussed by Mitra & Deshpande (1999). The phase angle, <j > , is varied between 
— 7 r and 7 r and divided into 500 bins with the flux contributions of core and conal components 
evaluated and summed. In order to see if the simulated pulsar is detected, the averaged flux, 
S ave , of the pulse profile is then compared to the flux threshold, S m i n , of each survey at 
its corresponding frequency. If detected, the pulsar is flagged as radio-loud otherwise it is 
radio-quiet. 

In the ACC study, pulsar surveys were all selected near a frequency of 400 MHz; hence 
they had no need to introduce any frequency dependence to the fluxes of the core and conal 
components. An important assumption made in ACC is that the ratio of the core to cone 
flux is given by r = yP -1 . In the present study, we have two groups of pulsar surveys. One 
group has frequencies near 400 MHz, while the other group has frequencies near 1400 MHz. 
Therefore, we have had to introduce the above frequency dependence to the spectra of the 
core and conal components with a ratio of the core to cone peak flux given by 


r = 10P' 1 


/ CX C ore d" 1 
\Q!cone d“ 1 


V 


,400 MHz, 


OLcore ~&cone. 
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( 8 ) 


where this ratio is similar to the ratio used in ACC at a constant frequency of 400 MHz. With 
this ratio, short period pulsars will have their radio fluxes dominated by the core component 
with a weak conal component depending on the viewing geometry. 

The luminosities of the core and cone components are expressed by 


L 


cone 


L 


core 


L 

(1 + r/r 0 ) ’ 

L 

(1 + r 0 /r)’ 


and 

where 
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^ cone ( &core + i \ /_ l / \ &core~&cone 

^ = Wone + l) V50 MHz/ ’ W 

where r/r 0 is the ratio of core to cone luminosities, and L is the total luminosity given by 

L = 3.4 x lO^P" 1 - 3 ^ 0 - 4 ( m jy • kpc 2 • MHz). (10) 

This luminosity is reduced by a factor of 60 from the one used in ACC, as discussed in 
section 5. Under this assumption, radio pulsars are believed to be standard candles with 
well-defined luminosities in terms of only the period and period derivative whose exponents 
in the expression above for L come from parameters in Table 1 of ACC for the first model 
assuming a braking index n = 3 with two velocity components. There is no dithering of the 
luminosity as in the case of Narayan & Ostriker (1990). Here the random viewing geometry 
accounts completely for the required dithering when the beam and viewing geometries are 
not included. However as discussed later in the text in section 5, in order to obtain a 
resonable birth rate and adequate agreement between the distributions of the distance, flux 
and dispersion measure, we have had to reduce the radio luminosity used in ACC by a 
substantial amount. 


3. Gamma-ray Emission Geometry 

For the geometry of the 7-ray beam, we have adapted the emission from the slot gap 
described in Muslimov & Harding (2003). The slot gap (Arons & Scharlemann 1979) is 
a narrow region between two conducting boundaries, the last open field line and the pair 
formation front, extending from the neutron star surface up to the light cylinder. Since 
the electric field is relatively small in the slot gap, primary particles accelerate more slowly 
and pair cascades form at altitudes of several stellar radii above the surface. We model the 
radiation from these pair cascades as having two components, curvature radiation of the 
primary electrons and synchrotron radiation from the electron-positron pairs. We obtain the 
number of curvature, Ncr, and synchrotron, N S r, photons emitted per primary particle by 
integrating the differential production rates over the available energy and over the pulsar 
phase angle <j>. Since we have not included relativistic effects such as aberration and time- 
of-flight delays in modeling the 7-ray beam, the caustic peaks as found by Dyks & Rudak 
(2003) and in outer gap models (e.g. Romani & Yadigaroglu 1995) do not appear in our 
calculations. 
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3.1. Curvature Radiation 

We assume that the curvature radiation takes place in the slot gap at the last open 
field line (see Figure 1 ) and is integrated from a lower 7 -ray threshold, E y , of 100 MeV to 
the curvature radiation critical energy, where E y and Ecr, are in m e c 2 units and having an 
angular distribution represented by the expression 


dN CR {> Ey, 0 7 ) 

d£t 


Ecr 

P 

r 


3e 2 


27rficsin0 

Hh ‘ 

m e cp 


(^) 1/3 [E&-BS?], where 


( 11 ) 
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97rr e 7; 


4cP 

r(l + cos 2 0) 3 / 2 

3 sin 0(1 + cos0) ’ 

CP . 2 n 
— sin 0 , 

27 r 




is the curvature photon critical energy and 


is the radius of curvature with 


where r is the radial distance on the last open dipole field line corresponding to magnetic 
colatitude, 0, and 0 7 = |0 is the photon emission angle, which is tangent to this field line. R 
is the neutron star radius taken to be 10 6 cm, r e is the Compton wavelength of an electron 
and 7 o is the initial Lorentz factor of the particle given by 


7. = for B^P-W < 1440, or (12) 

7 0 = 4 x 10 7 , for B\£ P~ 27 P > 3.0 in regime I 
7 0 = 1.4 x 10 7 P- 1/4 , for B^P- 27 ' 1 < 3.0 in regime II. 

where the first expression is for the case where the electric field is not screened by electron- 
positron pairs (Harding et al. 2002) and the second and third expressions are for cases where 
the electric field is screened by pairs in the unsaturated (I) and saturated (II) regimes (Zhang 
& Harding 2000). P is the period in seconds and P 12 is the magnetic field at the surface in 
units of 10 12 G. At a given inclination angle a, the line of sight angle, £, and phase angle, 
(j>, define a polar angle, 0 (through equation [3]), where the emission occurs tangent to the 
last open field line at a radial distance, r, from the stellar center, and the curvature emission 
rate per primary particle is given by equation ( 11 ) for dNcR/d£l. 
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3.2. Synchrotron Radiation 

The synchrotron radiation from cascade pairs takes place along the slot gap, beginning 
at an altitude Rmi n , where the first pairs are produced, and continuing out to R S r, the 
maximum radius at which pairs are produced (see Figure 1). R S r is determined as the 
altitude where the curvature radiation critical energy, Ecr, is equal to the photon escape 
energy, E esc (i.e. the minimum energy of photons capable of one-photon pair production). 
The pair-escape energy, E esc (r), in m e c 2 units, is given approximately by (see Zhang & 
Harding 2000 or Harding 2001) 


/ r p\V 2 / o 

E eS c(r) «518 ( — 1 max ( !, — I 


(13) 


The angles 0 m j n and 0 S are those corresponding to the radii Rmin and R S r along the 
field line at the edge of the slot gap, 


0s 



(14) 

(15) 


The corresponding photon emission angles are 0 7jm i n « §0 m in and 0 7iJ « § 0 S . The parameter, 
Rmin/R is set to 3.5 in all of our simulations, fixing the beginning of the emission zone to 
be 2.5 stellar radii above the surface (Muslimov Sz Harding 2003). We assume that the 
electron-positron pairs in the cascade have a spectrum N±( / y p ) = C p j~ p , extending from a 
minimum 'y min = E esc {R m m)/2 to a maximum at j max = E C r/ 2, where E esc (R min ) is the 
photon pair-escape energy at radius R m i n - 

The integral photon spectrum above energy E 7 of the synchrotron radiation from the electron- 
positron pairs with spectral index p, per primary particle, is 


N S r(> E y ) 


Esr 



where 




(16) 
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Esr is the critical synchrotron energy, in m e c 2 units, of pairs at their maximum energy 7 max 
and B' is the local field strength in units of the critical field 4.4 x 10 13 G. The pitch angle ip 
of the pairs is assumed to be that of the parent photon direction with respect to the local 
field at the pair production point, or 


simp = 
sin ip — 


wit’ for B ' < 01 ’ 01 

■J-, for B' > 0.1. 

Ecr 


(17) 


The spectral index of the pairs is given by p = 2a n - 1, where the spectral index of the 
photons, a n , is determined (Harding k Daugherty 1999, Zhang k Harding 2000) by the 
number of generations, n, of the pair cascade, 


Oi n 

n 


2 - 


2 — c*i 
2 n_1 ’ 


1 1 ^{Eesc/ E 0 ) 

ln(/c) 


(18) 


(19) 


where k = 3/64 and au = 5/3 for curvature radiation with losses. The cascade generation 
number is determined by E a , which is the curvature radiation critical energy at the initial 
colatitude angle 0 m j n of the radiation zone. 

We normalize the pair spectrum to the total cascade multiplicity 


N± ~ 


E 0 

20E esc ( < R m i n ) 


( 20 ) 


so that the normalization factor Cs for the synchrotron photon spectrum in equation (16) is 
given by 


C s = 


4tt 2 1-p (1 — p) 
3.3\/3 (3p — l) 


N± [Et P) 


EesciRru in) (1 " P) ]. 


( 21 ) 


The emission from the high altitude ( 2-4 stellar radii) cascades from the slot gap along 
the last open field line forms a broad, hollow-cone beam. The parameter representing the 
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longitudinal thickness of the slot gap is expressed, in units of the polar cap half-angle 0 min , 
as (Muslimov & Harding 2003) 


A£ = 4PP~ 4/7 for PP“ 2 4/7 < 0.075 or (22) 

Af = 0.3 for PB~ 2 4/7 > 0.075, 


where B x2 is in units of 10 12 G. The acceleration-cascade simulations indicate that the width 
of the slot widens as the pulsar ages and saturates at a value of approximately 0.3. As seen 
in Figure 1, the interior and exterior polar angles of the radiation from the slot gap at R min 
are described by the following expression 


aSG 



~ 0 7 ,min(l - A£), and 

min 


(23) 


We take the average opening angle of the cascade radiation from the slot gap between 
r = R m j n and r = Rsr as 8sg = (&m?n + 0~ U s)/2. We approximate the angular distribution of 
the synchrotron radiation component of the entire cascade between 6 min and 9 S as a hollow 
beam with a conal Gaussian of width equal to 


•u % 


= e SG - e 


iSG 

min 5 


(24) 


which is the full width at 1/e of the maximum. The integral photon spectrum above energy 
E 1 of the synchrotron radiation, per primary particle at a given polar angle is then given by 


dNsR^> 

dQ, 


N S r(> E 7 ) 

fe^ exp 


( 0 7 - 9sgY 

,.,2 


(25) 


The current of primary electrons in the slot gap that results in curvature and synchrotron 
radiation is limited to a fraction of the Goldreich- Julian current, n GJ , in the following m ann er 


h s iot gap = (1 - (1 - A 0 2 )n G J, where 


(26) 
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5.7 x 10 3l B 12 
n G j = pj “• 

( 27 ) 

This current multiplies the integral of he curvature and synchrotron emission per primary 
particle to give the total slot gap emission beam. 

The total flux due to both curvature and synchrotron radiation are calculated for a given 
phase angle, </>, which is related to 0 7 through equation (3), for a pulse profile with 500 bins 
of phase angle from — 7r to 7r. The average of the profile is obtained and compared to the 
appropriate instrumental flux threshold. If the average flux is above the threshold, the 7-ray 
pulsar is detected. This condition is tested independently of the radio flux and appropriate 
radio survey threshold, allowing us to designate the detected 7-ray pulsar as radio-loud or 
radio-quiet. 


4. Monte Carlo Simulations 

We discuss here some of the important changes that have been made to our Monte Carlo 
simulation code from our previous work in Gonthier et al. (2002). While we believe that 
is important to place the neutron stars at birth in spiral arms, we have not yet included 
the spiral arm structure into our simulations. As in Gonthier et al. (2002), we distributed 
pulsars at birth in the Galactic disk according to the prescription of Paczynski (1990). In 
a cylindrical coordinate system, the azimuthal angle, <f > , is randomly chosen between 0 and 
27r. The 2: distribution varies exponentially with distance from the plane, while the radial 
distribution peaks at 4.5 kpc and decreases exponentially from the center of the Galaxy. 
Given the initial position and velocity, the trajectory of each neutron star is evolved in the 
Galactic potential to the present. 


4.1. Comparison group of pulsars in the ATNF catalog 

In order to have a comparison group to normalize our simulation, we have selected 
pulsars from the Australian Telescope National Facility (ATNF) *. We chose pulsars within 
the Galaxy with periods greater than 30 ms and with positive period derivatives to obtain 
a comparison group of 978 pulsars detected by these nine surveys. Selecting pulsars with 


1 ATNF catalog is available at: http://www.atnf.csiro.au/research/pulsar/psrcat/. 
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periods greater than 30 ms insures that we have left out of our group most of the millisecond 
pulsars that have been recycled in binary systems. We are not currently simulating this class 
of pulsars since their evolution is more complicated. We have also not included the anomalous 
X-ray pulsars, the soft-gamma-ray repeaters, or pulsars in globular clusters in our comparison 
group. We run the Monte Carlo simulation until the code detects the same total number 
of radio pulsars as have been observed with the group of surveys. With this normalization 
a neutron star birth rate is predicted as well as the number of 7-ray pulsars detected by 
various instruments. However, in order to obtain smoother simulated distributions, we run 
the code for ten times the number of pulsars detected by the radio surveys and then normalize 
accordingly. 


4.2. Flux sensitivity of the Parkes multi-beam pulsar survey 

We have included the eight radio surveys described in Gonthier et al. 2002 along with 
the new PMBPS having an angular coverage of | b |< 5° and t — 260° to £ = 50° with an 
assumed geometric efficiency of 100%. In Gonthier et al. 2002, we calculated the m i nim um 
radio thresholds, S m i n , for the selected group of radio surveys using the Dewey et al. (1985) 
formula. We attempted to use the same formula for the PMBPS using the parameters in- 
dicated in Manchester et al. (2001). However, we found that an additional factor of ~ 2 
multiplying the limiting sensitivity is required to reproduce the S min curves in Figure 2 of 
Manchester et al. (2001) as shown here in Figure 2. The Dewey formula under predicts 
the S m in , and a more realistic treatment of narrow pulse widths (smaller duty cycles) in the 
Fourier search is not as optimistic as the Dewey formula (F. Crawford private communica- 
tion). As a result, we chose to evaluate the S min for the PMBPS using an IDL code (Crawford 
private communication) that we translated into C++ to incorporate into our Monte Carlo 
code that is called event-by-event. This routine was used to create Figure 2 in Manchester 
et al. (2001) and reproduced here in Figure 2, along with the S min curves predicted by the 
Dewey formula, with the extra factor of 2, for the indicated DMs. In our simulations, we 
have then scaled the limiting sensitivities as discussed in Manchester et al. (2001). 


4.3. Gamma-ray Thresholds 

We simulate the 7-ray pulsars detected by EGRET, AGILE and GLAST. If the sim- 
ulated 7-ray flux, obtained from the average flux in the pulse profile, is above a detector 
threshold, the pulsar is said to be a 7-ray pulsar detected by the corresponding instrument. 
We have included an all-sky sensitivity for both EGRET (I. Grenier, private communica- 
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tion) and AGILE (A. Pellizzoni, private communication) and are shown in Figure 3a and 
3b. For AGILE, we have three all-sky sensitivity maps representing the best, nominal and 
worst case scenarios. The one portrayed in Figure 3b is for the nominal case. For GLAST, 
we have used the following thresholds: in-plane (| b |< 10°) 5 x 1CT 9 photons/ (cm 2 • s), out- 
of-plane (| b |> 10°) 2 x 10 -9 photons/(cm 2 • s) (D. Thompson private communication), and 
for pulsed emission 5 x 10 -8 photons/(cm 2 • s) (S. Ritz private communication, McLaughlin 
& Cordes 2000). The above threshold for pulsed emission detection in a blind periodicity 
search is based on techniques used in periodicity searches of EGRET data (Mattox et al. 
1996, Chandler et al. 2001). 


4.4. New distance model 

We have incorporated the new electron density model of the Galaxy from Cordes & Lazio 
(2003) by calling the FORTRAN subroutines from within our code to calculate the dispersion 
measure (DM) of the simulated pulsar. The DM leads to a smearing of the pulse, affecting 
the flux threshold for radio detection. For comparison, we have recalculated the distance of 
the pulsars in the ATNF catalogue from the measured DM and the pulsars location using the 
new distance model. In Figure 4, we show the histogram for the logarithm of the absolute 
value of the difference between the distance obtained from the new distance model and the 
old distance model for our selected group of 978 pulsars. The distances obtained with the 
new distance model are about 20% smaller than those obtained with the old distance model. 
For pulsars in the catalogue whose best estimate distance is different than the one obtained 
using the old distance model, we have assumed that the distance was established by other 
methods and, therefore, is more reliable. 


4.5. Initial period distributions 

Recently various observations of young supernova remnants have been able to measure 
the speed of the expansion shell and the period and period derivative of the pulsar, thereby, 
determining the initial period of the pulsars. For example, X-ray Pulsars PSR J 1811-1925 
and PSR J0205+6559 (Gavriil et al. 2003) have been associated with the supernova remnants 
Gil. 2-0.3 and 3C 58, respectively, and may suggest that these pulsars where born with a 
period of ~ 65 ms. In contrast to our previous study of Gonthier et al. (2002) that used 
a constant birth period of 30 ms, we studied Gaussian and flat initial spin distributions 
to describe the initial period. We found that the overall population statistics are not very 
sensitive to the initial spin distribution and only affects the short-period population of pulsars 
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in the P — P diagram. While significant progress is being made in deducing the initial period 
of pulsars, the shape of the distribution is not well defined at the present. We have concluded 
that a flat distribution from 0 to 150 ms accommodates the observations and have used this 
distribution in this study. 


4.6. Decay of the Magnetic Field 

We continue to be steered in the direction of incorporating the decay of the magnetic 
field in order to achieve better comparisons. Originally, we included eight radio surveys in 
Gonthier et al. (2002) with 445 detected radio pulsars to compare to our simulated results. 
We used one Gaussian to describe the primary magnetic field distribution with a single decay 
constant. The PMBPS (Manchester et al. 2001) has discovered many more pulsars, many 
of which are young, distant pulsars with high radio luminosities. The current pulsar catalog 
now has 1412 radio pulsars (http://www.atnf.csiro.au/people/pulsar/catalogue). With the 
PMBPS, we have a selected group of 978 detected radio pulsars with many more high field 
pulsars. A single Gaussian would result in too many low field pulsars. In order to simulate 
these high field pulsars, we found it necessary to use two Gaussian distributions to skew the 
distribution towards high fields. 

The pulsars surface magnetic field distribution at birth is represented by the sum of two 
log-normal Gaussian distributions expressed as 

Pb = Yj Aie ' (1 ° 8 B_1 ° g Bi? ‘ a ' ’ (28) 

2=1 

where the parameters are indicated in Table 1. 


Table 1 



Parameters for Initial 

Magnetic Field Distribution 

i 

Ai log Bi 

<7» 

i 

0.6 12.75 

0.4 

2 

0.3 13.0 

0.7 


While there are two Gaussians describing the initial field distribution of the pulsars at birth, 
the second Gaussian with a higher mean field merely skews the distribution towards higher 
magnetic fields and does not necessarily suggest two groups of pulsars with different field 
characteristics. Using a single broader Gaussian would result in too many lower field pulsars. 
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The birth rate is assumed to be constant during the history of the Galaxy (at least back 
to 10 9 years in the past); therefore, we randomly select the age of the pulsar from the present 
to 10 9 years in the past. We assume a dipole spin down with a decaying magnetic field having 
a time constant, t d , following Gonthier et al. (2002). In Figure 5, we present the period 
derivative versus the period of simulated pulsars with the indicated time constants from 10 8 
to 5 x 10 5 years. Indicated in the figure are lines of constant field (calculated according to 
Shapiro & Teukolsky, 1983) and pulsar age assuming dipole spin-down of a constant field as 
well as the curvature radiation (CR) and nonresonant inverse Compton scattering (NRICS) 
death lines (Harding et al. 2002). Without field decay or a large decay constant like 10 8 
years, pulsars will move from their short periods to longer periods along constant or nearly 
constant field lines and pile up near the NRICS line. 

In order to populate, without field decay, the diagram in a region of small period deriva- 
tives (5 x 10 -18 ) and medium periods (0.5 s), many more short period pulsars would populate 
the lower left region of the diagram where none are observed (in Figure 8). Decreasing the 
decay constant, produces the upside-down pear shaped distribution seen in the distribution 
of detected pulsars and populates the high field region above 5 x 10 13 G. Unless one can 
alter the period and period derivative dependence of the radio luminosity significantly in 
a manner with more than just a simple power law, we find that field decay is required to 
reproduce the distribution. In subsequent simulations, we have adopted a value of 2.8 Myr 
for the decay constant. 


4.7. Supernova kick velocity distribution 

A number of studies disagree on the initial 3-D velocity distribution of neutron stars 
at their birth, possibly the result of an asymmetric supernova kick, typically described by a 
Maxwellian distribution. Lorimer, Bailes & Harrison (1997) obtained a velocity distribution 
with a mean velocity of ~480 km/s similar to a previous study obtaining a mean of ~450 km/s 
(Lyne & Lorimer 1994), yet significantly larger than most previous studies of pulsar statistics 
that required space velocities of ~150 km/s. Hansen & Phinney (1997) concluded that a 
mean velocity ~250-300 km/s best described their study. Though Hartman et al. (1997) 
did not use a Maxwellian distribution, they obtained a distribution with a mean velocity of 
380 km/s. Gonthier et al. (2002) also did not use a Maxwellian distribution and found a 
distribution with a mean velocity of 170 km/s. It is clear that the velocity distribution that 
one obtains depends heavily on the many other assumptions that go into the model, such as 
the radio luminosity and radio beam geometry. The brighter the radio pulsars, the greater 
the distance at which they are detected, resulting in a broader distribution of distance, z, 
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from the plane of the Galaxy, requiring a smaller mean velocity to improve the agreement 
with 2 distribution of detected pulsars. 

In this study we have adopted the luminosity model of ACC, and so we must also adopt 
their kick velocity model. We chose to follow their two-component velocity distribution, 
which is Maxwellian in velocity with characteristic widths of 90 and 500 km/s and given 
by the equation (1) in ACC. In the ACC model, the two-component velocity model was 
preferred over the single component model. The fraction of the neutron stars with a width 
of 90 km/s is 40% and with a width of 500 km/s is 60%, leading to an average velocity of 
rsj 540 km/s. 

In Figure 6, we show the 2 distribution above the Galactic disk for the detected pulsars 
(shaded histogram) and for the simulated pulsars (regular histogram). Under the assump- 
tions of the model, the predicted distribution is a little wider that the one for detected pulsars, 
having scale heights of 152 pc and 182 pc, respectively. We realize that many assumptions 
in our model are interrelated and decreasing the overall radio luminosity of the ACC model 
leads to a difference in the z distribution as pulsars are radio dimmer and must be closer to 
be detected. We have chosen to keep the velocity model and overall radio luminosity model 
of ACC making as few necessary changes to the ACC model as needed. 


4.8. Reduction of the radio luminosity 

Using the radio luminosity of ACC, we find that the simulated radio pulsars are too 
bright, with too many distant pulsars being detected and predicting a neutron star birth 
rate of 0.11 per century, with no 7-ray pulsars predicted to be detected by EGRET. In 
ACC all the pulsar surveys used in the study were at frequency near 400 MHz. In the set 
of surveys chosen for this study there are two groups with one having frequencies near 400 
MHz and the other group with frequencies near 1400 MHz. Since the PMBPS S m { n is best 
accounted for in our simulation and this survey at 1374 MHz detected most of the pulsars 
observed in the Jodrell Bank 2 survey at 1400 MHz and the Parkes 1 survey at 1520 MHz, 
we selected only the PMBPS pulsars to represent the high frequency (HF) surveys, while the 
other selected surveys in our study near 400 MHz represent the low frequency (LF) surveys. 
Focusing primarily on the distributions of the pulsar distance, DM and flux at 400 MHz and 
1400MHz, we first set the spectral indices (preserving 0.5 between core and cone indices) to 
give the same birth rate for these two frequency groups, then we set the over all luminosity 
to given a reasonable birth rate of ~ 1.5 neutron stars per century. We obtain the same birth 
rate for each group with spectral indices of a core = —2.1 and a cone = —1.6. These spectral 
indices describe the primary spectra of the radio pulsars before the selection effects of the 
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characteristics of the chosen radio surveys. 

In our simulation, we calculate the fluxes at the frequency of each of the surveys in our 
selected group as well as at 400 and 1400 MHz by averaging the pulse profile for a given 
random viewing geometry for each pulsars. From the calculated S400 and S1400, we obtain 
a spectral index for each pulsar. In the simulation with all the surveys, we simulated 9780 
(ten times the number in our select group of surveys to improve statistics) radio pulsars as 
detected by these surveys and find an average spectral index of -1.8 and a standard deviation 
of 0.2. Lorimer et al. (1995) measured spectral indices of 280 pulsars by measuring fluxes 
at radio frequencies between 408 and 1606 MHz with the distribution having a dependence 
on the characteristic age of the pulsar with a dependence, a = — 1.7 + 0.2 log with 

a standard deviation of 0.6 with respect to this dependence. Maron et al. (2000) extended 
to lower and higher frequencies the study of Lorimer et al. (1995) and obtained an average 
spectral index of -1.8 with a standard deviation of 0.2. 

To show the overall effect of the reduction of the radio luminosity used in ACC, we show 
in Figure 7 the comparisons of the distributions for the pulsar distance, flux at 1400 MHz and 
DM for 620 pulsars detected with only the PMBPS and an equal number simulated pulsars 
for this survey alone. Given that we are calculating the S m i n for the PMBPS according to the 
formulation used in Manchester et al. (2001), we believe that the S m i n for this survey is the 
best described in the group of our selected surveys, and with all pulsars detected at ~1400 
MHz sky temperature effects are minimized, reducing further uncertainties. Indicated in 
Figure 7 are the factors, f Te d, used to reduce the radio luminosity of the ACC model and the 
resulting distributions. Since the ACC model studied pulsar surveys at 400 MHz, the factor 
fred represents a reduction of the ACC 400 MHz luminosity by this factor. As the radio 
luminosity is reduced, the comparison of the distances and DMs improves, but disagreement 
increases between the distribution of the simulated and detected radio flux at 1400 MHz. 
There are significantly more pulsars simulated with lower radio fluxes than those detected 
by this survey, suggesting that perhaps certain aspects of the emission geometry are still not 
adequately described. We find that the shape of the flux distribution is not very sensitive to 
the features of the conal geometry, such as the radius and width. We also wanted to predict 
a reasonable neutron star birth rate per century, which varies from 0.6 at f Te d = 20 to 1.6 
at f rei i = 80. We chose to reduce the radio luminosity by a factor of f Te d — 60 in subsequent 
simulations in the following figures compromising between good agreement of the distances 
and DMs and less desirable agreement of the fluxes. These factors predict a birth rate of 
1.38 neutron stars per century for the case of all nine radio surveys. 
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5. Results 

In Figure 8, we present the distribution of our select group of 978 detected pulsars in 
the ATNF pulsar catalog (8a) and the same number of simulated pulsars (8b) as a function 
of period derivative and period. The solid lines represent constant dipole magnetic field 
of 10 11-14 G, and the dashed lines correspond to the curvature and nonresonant inverse 
Compton scattering death lines. As dotted curves, we show paths in the diagram for 4 
pulsars assuming a field-decay constant of 2.8 Myr, all with ages of 10 7 years, and with the 
indicated initial magnetic fields in B X2 units. We show the age lines assuming field decay as 
dot-dashed lines. Due to field decay, these lines are very different from the characteristic age 
lines, with the oldest pulsars being 10 Myr rather than 1 Gyr assuming no field decay. As 
indicated in Table 1, we chose to represent the initial magnetic field distribution with two 
Gaussians. The initial magnetic fields of the four evolutionary paths portrayed in Figure 8 
were chosen to represent the means of the two Gaussians, the higher mean plus the high 
field component’s width and the lower mean minus the lower field component’s width. As a 
result, most of the simulated pulsars will lie within the lowest and highest paths. Due to a 
2.8 Myr decay constant, the knee-like, curved portions in these paths begins before 1 Myr 
and become vertical after a few Myr. The paucity of pulsars in the two regions indicated 
with the circle and LB and HB for low and high field (Figure 8a) can be explained in terms 
of decay of the surface magnetic field on a time scale of 2.8 Myr. While it may be true that 
one can tailor a radio luminosity law to account for the observed distribution, field decay 
leads naturally to the funnel-shaped distribution of the detected pulsars. 

In Figure 9, we present the comparisons between detected and simulated statistics for 
the indicated parameters. The shaded histograms represent the distributions of the detected 
pulsars while the regular histograms correspond to the model simulations. The model sim- 
ulation over-predicts the number of pulsars with short periods, large period derivatives and 
larger distances. Improved comparisons can be obtained by making the radio luminosity 
not as strongly dependent on the period by decreasing the exponent from -1.3 to -0.9 and 
that of the period derivative from 0.4 to 0.3. The main discrepancy lies in the comparisons 
of the distance distributions and flux distributions (Figure 5). The agreement in the flux 
distributions for only the PMBPS, shown back in Figure 5, is better with a reduction factor 
of fred — 20 making the pulsars more radio bright. However, the agreement between the dis- 
tance distributions is significantly worse and the birth rate is too low. Perhaps the inability 
to find agreement in the radio flux and distance distribution may stem from our underlying 
assumption that radio pulsars are standard candles or may be indicating that spiral arm 
structure may have to be used for the birth location of neutron stars in order to improve the 
agreement between the simulated and observed distributions. 
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In Figure 10, we show the distributions in Galactic coordinates as Aitoff projections for 
978 pulsars detected (10a) and simulated (10b). The strong contribution to the Galactic 
disk is due primarily to the PMBPS, adding nearly half of the total number of pulsars that 
are unique to this survey. 

We also simulate the 7-ray pulsar detections by EGRET, AGILE and GLAST, using 
the assumptions discussed earlier, and independently simulate detection of radio pulsars. 
However, the 7-ray pulsars are flagged as radio-loud if their fluxes are higher than the min- 
imum sensitivities of the select group of surveys; otherwise they are radio-quiet. Therefore, 
we can predict the number of radio-loud and radio-quiet 7-ray pulsars detected as point 
sources by each of the three instruments. For known radio pulsars, with measured periods 
and period derivatives, the 7-ray instrument can detect them as point sources and obtain a 
pulsed detection through reliable epoch folding. GLAST will have the sensitivity and ability 
to perform blind period searches and detect pulsation without a radio ephemeris. Table 2 
indicates the simulated pulsar statistics for the radio-quiet and radio-loud 7-ray pulsars as 
detected by the three instruments. 


Table 2 


Simulated Pulsar Statistics 



All Nine Surveys 
(978 radio pulsars) 

All Surveys Excluding PMBPS 
(546 radio pulsars) 

Birth Rate 

1.38 

1.46 


Radio Quiet 

radio-loud 

radio-quiet 

radio-loud 

EGRET 

7 

19 

10 

15 

AGILE 

13 

37 

19 

36 

GLAST 

276 

344 

436 

209 


To improve the statistics, we run the simulation until the number of radio pulsars 
detected by the chosen set of radio surveys is equal to ten times the number of actual detected 
radio pulsars by those surveys. We then renormalize by dividing by a factor of ten. EGRET 
observed all of the radio pulsars that were detected by all the radio surveys excluding the 
PMBPS, out of which EGRET detected 8 7-ray pulsars including Vela, the Crab, B1951+32, 
B1706-44, B1055-52, B0656+14, J1048-5832, and the radio-quiet 7-ray pulsar, Geminga plus 
a couple of other candidate pulsars bring the total to perhaps 12 (Kanbach 2003). Excluding 
the PMBPS, our simulation predicts that EGRET should have seen 15 radio-loud and 10 
radio-quiet 7-ray pulsars and a neutron star birth rate of 1.46 per century. With all nine 
surveys, the predicted birth rate of 1.36 is 5% smaller and, therefore, the total number of 7- 
ray pulsars has also dropped by 5% (GLAST) within statistics. The EGRET Third catalogue 
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(Hartman et al. 1999) contains 170 unidentified point sources, some of which are expected to 
be pulsars. Sensitive searches performed with Chandra and the Parkes multi-beam telescopes 
have resulted in a few new pulsars within the error boxes of the unidentified EGRET sources 
(Halpern et al 2001; D’Amico et al. 2001). Correlating the positions of the radio pulsars 
detected in the PMBPS with the EGRET unidentified sources, Torres et al. (2001) found 
14 positional coincidences. With the nearly completed PMBPS, Kramer et al. (2003) found 
about 38 positional coincidences, and they determined that 19 ± 6 are statistically likely to 
be real associations. So it would seem then that adding the PMBPS should convert radio- 
quiet 7-ray pulsars into radio-loud 7-ray pulsars as detected by EGRET. This is clearly the 
case for GLAST and there is a significant conversion of radio-quiet to radio-loud pulsars 
for AGILE and EGRET when the PMBPS is added. Since our present simulation over- 
predicts the number of short period, energetic pulsars, mostly detected by the eight surveys 
(without PMBPS) prior to EGRET, our results are over-predicting the number of radio- 
loud pulsars detected by EGRET from these surveys. However, the number of radio-loud 
pulsars predicted for all surveys including PMBPS is consistent with the number of plausible 
coincidences found by Kramer et al. (2003). 

In Figure 11, we present the positions in the P - P diagram of the known 7-ray pulsars 
detected by EGRET (10a) and of those simulated for EGRET (10b), AGILE (10c) and 
GLAST (lOd) where radio-loud 7-ray pulsars are shown as solid circles and radio-quiet 
7-ray pulsars are shown as crosses. Younger pulsars have higher 7-ray luminosities that 
decrease as they approach the curvature death line where curvature radiation 7-rays can 
no longer produce electron-positron pairs. The 7-ray luminosity decreases significantly for 
older pulsars below the curvature radiation death line (Harding et al. 2002), where the 
main mechanism of pair production is via inverse Compton scattering of the thermal soft X- 
rays from the stellar surface. Pulsars below the nonresonant Compton scattering death line 
are unable to produce pairs and become radio-quiet. Our simulations predict that GLAST 
will detect 276 radio-quiet and 344 radio-loud 7-ray pulsars; a significant improvement over 
EGRET. With GLAST ability to perform blind period searches, we predict that out of 
the 276 radio-quiet 7-ray pulsars, 17 pulsars will be detected as pulsed sources. AGILE, 
scheduled to launch before GLAST, should detect 13 radio-quiet and 37 radio-loud 7-ray 
pulsars using the nominal sensitivity map. The best and worse case-maps predict 24 and 9 
radio-quiet and 64 and 27 radio-loud 7-ray pulsars, respectively. As all-sky sensitivity and 
threshold maps are improved, these numbers will vary somewhat. 

We present in Figure 12 the Aitoff projections of the 7-ray pulsars observed by EGRET 
and those simulated for EGRET, AGILE and GLAST. Most of the pulsars detected by 
GLAST will be young pulsars with ages ~ 10 5 years and within 500 pc of the Galactic disk. 
The asymmetric distribution of radio-loud and radio-quiet 7-ray pulsars is a result of the 
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PMBPS detecting radio pulsars from i = 260° to l = 50°, mostly on the right side of the 
figure. 

In order to further explore the effect of the parameterization of the beam geometries of 
the radio and 7-ray emission, one has to study the pulse profiles and the correlation between 
the radio profile and the 7-ray profile. In Figure 13, we present a select pair of examples 
of simulated radio and 7-ray profiles for two radio-loud pulsars “detected” by EGRET with 
periods of 528 ms (13a) and 41 ms (13b). The two pulsars also have similar impact angles of 
4.1° (13a) and -5.5° (13b). The radio profiles presented here are for 400 MHz and the 7-ray 
profiles are for greater than 100 MeV. The ACC model predicts ratios of the radio core-to- 
cone flux of 13 (Figure 13a) and 163 (Figure 13b), suggesting that the conal contribution to 
the flux is minor. However, the viewing geometry is important also. The pulse profile in the 
line-of-sight intersects the radio cone and the outer portion of the radio core, displaying two 
radio peaks, but only one curvature radiation 7-ray peak. Due to the large impact angle, 
the radio core does not significantly contribute to the overall profile as this is a fairly long 
period pulsar. The other pulsar in Figure 13b has a similar impact parameter but with a 
shorter period, therefore the radio core clearly dominates the profile and the two conal peaks 
from the 7-ray pulse are manifested. Since in this model the 7-ray emission originates within 
2.5 stellar radii of the surface and the radio core emission is believed to come from similar 
altitudes above the stellar surface, there would be little aberration or time delay between 
them and they should appear correlated with the radio core peak being in phase with the 
single 7-ray peak, as in Figure 13a, or being in between the two 7-ray peaks, as in Figure 
13b. On the other hand the radio cone beam is believed to arise from a higher altitude region 
and the effects of aberration and time delay, not included in our model, might shift the peaks 
to earlier phase relative to the 7-ray or radio core in the pulse profile. This effect is actually 
observed in a number of triple radio profiles, where the core component lags the center of 
the cone (Gangadhara & Gupta 2001). In the ACC model with core dominated short-period 
pulsars, we find that 54% of the EGRET radio-loud 7-ray pulsars have two 7- ray peaks in 
the pulse profile and are all core dominated, exhibiting a single (core) radio peak in the 
profile. The other 46% have a single 7-ray peak with a variety of radio profiles, including 
a single core peak, core and conal peaks, and two conal peaks. However, the fact that all 
two-peaked 7-ray profiles are core dominated is contrary to the few EGRET 7-ray pulsars, 
typically consisting of two 7-ray peaks with a single radio peak leading the 7-ray peaks 
(Thompson et al. 1997). This significant discrepancy between the observed and simulated 
radio and 7-ray profile correlations either questions the parameterization of the core to cone 
radio flux in the ACC model suggesting that these short period pulsars are different, having 
their profiles dominated by cone emission, or questions the polar cap 7-ray emission model 
we have used. 
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By contrast, several recently discovered, young X-ray pulsars show single broad peaks 
in their profiles and also single radio peaks in phase with the X-ray peaks (Lorimer 2003 
and references therein). They also have very low radio flux and luminosity. Harding et 
al. (2003) have shown that the radio characteristics can be understood in the ACC model 
if we are viewing the edge of the cone beams at large impact parameter, so that the cone 
beam appears single and the core component is not seen. In the polar cap model emission 
geometry, the high-energy emission peak would be in phase with the single radio peak, as 
observed. A number of these are coincident with EGRET sources and AGILE may measure 
the 7-ray profiles. 


6. Discussion 

We have included radio and 7-ray beam geometries into our Monte Carlo code that 
simulates population statistics for radio and 7-ray pulsars, making predictions for the number 
of radio-quiet and radio-loud 7-ray pulsars detected by the instruments: EGRET, AGILE 
and GLAST. The radio beam geometry is tailored after the phenomenological model of ACC 
with slight modifications to include the radius-to-frequency mapping of Mitra & Deshpande 
(1999). In the ACC model, radio pulsars are assumed to be standard candles with their 
radio luminosity described by a simple power law in period and period derivative. The 7-ray 
beam geometry has been derived from the theoretical work of Muslimov &: Harding (2003) 
describing the emission from the slot gap in the polar cap model. These enhancements 
are significant improvements over our previous studies of Gonthier et al. (2002). We have 
added the PMBPS to our select group of radio surveys, we have used the new distance 
model of Cordes & Lazio (2003), we have added all-sky threshold maps for EGRET (Grenier 
private communication) and AGILE (Pellizzoni private communication), and we have used 
realistic S m i n thresholds for the PMBPS (Crawford private communication). Neutron stars 
are assumed to be born with a constant birth rate, with a Galactic distribution described by 
Paczynski (1990) with a flat distributions in their initial period and a Gaussian distribution 
in their magnetic field, and their trajectories are evolved to the present within the Galactic 
potential (Paczynski 1990). We find better agreement when the magnetic field is assumed to 
decay exponentially with a decay constant of 2.8 Myr, which is somewhat shorter than the 
decay constant in Gonthier et al. (2002) that assumed an entirely different beam geometry 
and radio and 7-ray luminosity models. 

In order to obtain agreement between observed and simulated radio pulsar distance 
distributions, a reasonable birth rate of 1.45 neutron stars per century and a reasonable 
number of 7-ray pulsars observed by EGRET, we find it necessary to significantly reduce 
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by a factor 60 the overall radio luminosity of ACC model. The major problem with our 
radio emission model is that we are not able to simultaneously fit the distance distribution 
and the radio flux distribution as seen in Figure 5. We found that slight adjustments to the 
radio beam geometry do not affect these distributions. This inadequacy perhaps challenges 
some major assumptions of the model such as that radio pulsars are standard candles or 
that short period pulsars are core dominated. There is a lack of theoretical insight into the 
mechanism for the radio luminosity. Randomizing the radio luminosity about the expected 
value does indeed achieve some improvement in the flux distribution. The disagreement 
between the inferred and the simulated distance distributions may be a result of assuming a 
smooth distribution of neutron stars in the Galactic plane, rather than taking into account 
the spiral structure of the Galaxy. It is very apparent that, most pulsars are detected in 
spiral arm regions. 

Our simulation is normalized to the number of radio pulsars observed by the selected 
group of surveys and, therefore, a neutron star birth rate is predicted, as well as the number 
of 7-ray pulsars observed by each instrument as shown in Table 2. We expected by adding the 
PMBPS that the ratio of radio-quiet 7-ray pulsars to radio-loud 7-ray pulsars would decrease 
for EGRET given the positional coincidences found between the Parkes multi-beam radio 
pulsars within EGRET error boxes. Recently Kramer et al. (2003) estimated about 19 ± 6 
associations, which suggests about 20 radio-quiet 7-ray pulsars detected by EGRET that 
become radio-loud when the PMBPS is included. The simulated ratio of radio-loud to radio- 
quiet 7-ray pulsars decreases by 45% from 2.7 (with PMBPS) to 1.5 (without PMBPS) for 
EGRET. For AGILE there is also a similar decrease of 34% from 2.9 (with PMBPS) to 1.9 
(without PMBPS), but very evident for GLAST as the ratio decreases significantly by 62% 
from 1.25 (with PMBPS) to 0.5 (without PMBPS). 

Our model assumes that pulsars are evenly distributed azimuthally throughout the 
Galactic disk peaking radially at 4.5 kpc and falling off exponentially with distance from the 
center. Our model does not assume the spiral structure of the Galaxy where the pulsars are 
indeed born. The location of most of the detected pulsars is clearly correlated with the spiral 
arms of the Milky Way. Perhaps the difficultly of our model in reproducing the distance and 
flux distributions is a result of the lack of including this strong correlation of the location of 
pulsars and the spiral arms. 

We believe that the fact that we cannot reproduce the decrease in the ratio of radio- 
loud to radio-quiet 7-ray pulsars for EGRET points to a significant limitation of our overall 
model. The fact that we are not able to simultaneously account for the detected radio flux 
and inferred distance distributions and that the radio pulse profiles do not correlate with 
the 7-ray profiles as those detected by the EGRET, also suggest that we are not adequately 
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accounting for everything. The problem with the pulse profiles perhaps implies that the 
parameterization of the ratio of radio core to cone with period in the ACC model does not 
apply for 7-ray pulsars. Recently Crawford, Manchester & Kaspi (2001) reported that six 
young pulsars have large linearly polarized profiles characteristic of conal emission. Yet in 
the ACC model, young pulsars with short periods would be characterized by significant core 
emission. In order to avoid the core component, the impact angle must be fairly large to 
observe a single peak from the edge of a conal beam. The pulsar J1105-6107 has a period 
of 0.063 s, and the ACC model predicts an integrated flux core to cone ratio of about 60 
making it difficult to avoid seeing the core beam when a significant mean conal flux of 1.0 
mJy is detected. Perhaps partial conal emission becomes dominant for short period pulsars 
as suggested by Manchester (1996). 

The correlation between the radio and 7-ray beam profiles and the shapes of the profiles 
are sensitive to the geometry of the beams as well as the viewing geometry. In the ACC 
model that we have used in our simulations, short period pulsars have fluxes that are core 
dominated. All two-peaked 7-ray pulse profiles display a single radio core peak in the radio 
pulse profile. In order to see a cone dominated radio profile, the impact parameter must be so 
large that single-peaked 7-ray profiles are seen. This appears contrary to the features of the 
most of the 7-ray profiles detected EGRET and again raises the issue of whether short period 
pulsars are core dominated and perhaps are really partial cone dominated as suggested by 
Manchester (1996) and, more recently, by Crawford, Manchester & Kaspi (2001). 

We conclude that a better radio beam model is required to account for the observed 
characteristics of radio pulsars. In addition, recent theoretical work by Dyks & Rudak 
(2003) and Muslimov &: Harding (2004) indicates the importance of including the caustic 
component of 7-rays in the emission geometry. We hope that in the near future, we will be 
able to include more realistic emission geometries of the radio and 7-ray beams. 
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Fig. 1. — Geometry of the slot gap emission relative to the magnetic moment vector, fi. The 
pair formation front (PFF) curves upward near the last open field line and forms the inner 
boundary of the slot gap with width A£ (in units of polar cap half- angle). R min and Rsr 
are the minimum and maximum radii of pair synchrotron radiation. 0^*? and 6^. are the 
tangent angles to field lines at the inner and outer edge of the slot gap at R m in . 
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Fig. 2. — The minimum radio flux, S min , detectable by the Parkes multi-beam pulsar survey 
for the indicated DM (in pc cm -3 ) as a function of period, calculated as in Manchester et 
al. (2001) (solid curves) and calculated using the formula in Dewey et al. (1985) (dotted 
curves) with an additional factor of 2 (see text) . 
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Fig. 3. — Gamma-ray flux threshold maps in Galactic latitude versus longitude for EGRET 
(I. Grenier, private communication) in (a) and for AGILE (A. Pellizzoni, private communi- 
cation) in (b) for the nominal case. Gray scale is in units of photons s -1 cm -2 . 
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Fig. 5. — Distributions of simulated radio pulsars as a function of the period derivative and 
period for the indicated decay constants of the magnetic field. Dashed lines represent the 
death lines for curvature radiation (upper) and for nonresonant inverse Compton scattering 
(lower). Dot-dashed lines represent the indicated pulsar age (in yr) and solid lines represent 
the indicated magnetic surface field strength (in units of G) assuming a constant dipole 
spin-down field. 
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Fig. 7. — Distributions of the distance, flux at 1400 MHz and dispersion measure for detected 
pulsars in the PMBPS (solid histograms) and simulated pulsars (regular histograms) for the 
indicated factor, / re <j, reducing the radio luminosity of the ACC model. 
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Fig. 8. — Distributions of (a) detected and (b) simulated pulsars with a field-decay constant 
of 2.8 Myr. Dashed lines represent the death lines for curvature radiation (upper) and for 
nonresonant inverse Compton scattering (lower). Dot-dashed lines represent the indicated 
pulsar age assuming a field decay of 2.8 Myr and solid lines represent the indicated magnetic 
surface field strength assuming a constant dipole spin-down field. Dotted curves for the 
indicated field strengths assume a decay constant of 2.8 Myr of the magnetic field. Two 
regions of pulsar paucity are indicated by circles - one at low field (LF) and another at high 
field (HF). 
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Fig. 9. — Distributions of various pulsar characteristics indicated as shaded histograms 
(detected pulsars) and plain histograms (simulated pulsars). 



Fig. 10. — Aitoff plots of the (a) detected and (b) simulated radio pulsars. 
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Fig. 11 . — Distributions of radio-quiet (crosses) and radio-loud (solid dots) 7 -ray pulsars (a) 
detected by EGRET and simulated for (b) EGRET (c) AGILE and (d) GLAST, assuming 
a field-decay constant of 2.8 Myr. Dashed lines represent the death lines for curvature 
radiation (CR) and for nonresonant inverse Compton scattering (NRICS). Dot-dashed lines 
represent the indicated pulsar age assuming a field decay of 2.8 Myr and solid lines represent 
the indicated magnetic surface field strength assuming a constant dipole spin-down field. 
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Fig. 12. — Aitoff plots of radio-quiet (crosses) and radio-loud (solid dots) 7-ray pulsars (a) 
detected by EGRET and simulated for (b) EGRET (c) AGILE and (d) GLAST, assuming 
a field-decay constant of 2.8 Myr. 






